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Abstract 

We determine an asymptotic expression of the blow-up time tcoii for self-gravitating 
Brownian particles or bacterial populations (chemotaxis) close to the critical point. We 
show that tcoii = t*{ri — r]c)~^^'^ with t^, = 0.91767702..., where rj represents the inverse 
temperature (for Brownian particles) or the mass (for bacterial colonies), and rjc is the 
critical value of r] above which the system blows up. This result is in perfect agreement 
with the numerical solution of the Smoluchowski-Poisson system. We also determine the 
asymptotic expression of the relaxation time close but above the critical temperature and 
derive a large time asymptotic expansion for the density profile exactly at the critical 
point. 

1 Introduction 

Recently, several papers have focused on the blow-up properties of a cluster of self-attracting 
particles |Ill2lElllllSlEllIllHlinilinilIIllI21- This concerns in particular the "isothermal 
collapse" JSj of self-gravitating Brownian particles and the "chemotactic aggregation" jT3] of 
bacterial populations in biology. These systems have a similar mathematical structure and 
they are described, in a first approximation, by the Smoluchowski-Poisson system jTHj. The 
Smoluchowski equation is a particular Fokker-Planck equation in physical space involving a 
diffusion due to Brownian motion and a drift. In the standard Brownian theory developed by 
Einstein, the drift is due to an external potential ^Hl- Alternatively, we can consider the more 
complicated situation in which the potential is produced by the particles themselves. If we 
assume a Newtonian type of interaction, we have to couple the Smoluchowski equation to the 
Poisson equation. 

It can be shown that the Smoluchowski-Poisson system decreases a Lyapunov functional 
F[p\ which can be interpreted as a Boltzmann free energy J7j. Furthermore, the type of 
evolution depends on the value of a dimensionless parameter 77. For r] < r]c = 2.51755132..., the 
Smoluchowski-Poisson system evolves to an equilibrium state which minimizes the free energy 
at fixed mass. It corresponds to an isothermal distribution of particles similar to isothermal 
stars and isothermal stellar systems flHJ- On the contrary, for rj > rjc, there is no possible 



equilibrium state and the system blows up. This is the case for self-gravitating Brownian 
particles below a critical temperature Tc and for bacterial populations above a critical mass 
Mc. It is found that the collapse is self-similar and that it develops a finite time singularity, i.e. 
the central density becomes infinite in a finite time tcoii [3 El- Then, a Dirac peak is formed in 
the post-collapse regime [TT] . 

It is clear that the collapse time tcou depends on the distance to the critical point rjc and 
that it should diverge as t] ^ rjc. More precisely, by using heuristic arguments, it is argued 
in [7] that tcou ~ t*(?7c — 77)"^^^. This scaling is consistent with numerical simulations of the 
Smoluchowski-Poisson system. However, the approach of |7j is qualitative and does not provide 
the numerical value of t^. 

After some introductory material presented in Section 2, we develop a systematic procedure 
in Section 3, inspired from |T^, which confirms the scaling law tcoii ~ t^(ric — r])~^^'^ and leads to 
the explicit value of t^,. In section 4, we study the relaxation time r to equilibrium below rjc and 
show that it diverges like r ~ Crj{T]—r]c)~^^'^, where Cr, is given explicitly. In Section 5, we derive a 
systematic large time expansion of the density profile exactly at rjc. In particular, we show that 
the central density approaches its equilibrium value according to po — pij = 0,t) ~ where 
Cp is a universal constant which is explicitly computed. Finally, in Section 6, we determine the 
pulsation period of bounded isothermal spheres described by the Euler- Jeans equations close 
to the critical point. 



2 Self-gravitating Brownian particles and bacterial pop- 
ulations 

2.1 The model equations 

In the overdamped regime, a system of self-gravitating Brownian particles is described by the 
N coupled stochastic equations 

(1) ^ = -/iV.f/(ri,...,r^) + V2D Ri(t), 

where p = 1/.^ is the mobility is the friction coefficient), D is the diffusion coefficient and Ri(t) 
is a white noise satisfying (Rj(t)) = and {Ra,i{t)Rbj{t')) = 6ij6abS{t — t'), where a,b = 1, 2, 3 
refer to the coordinates of space and i,j = 1, to the particles. We define the temperature 
T = 1/13 through the Einstein relation p = D(3. The particles interact via the potential 
t/(ri, Tjv) = Ylii<j'^i^i^^j)- In this paper, u{Yi — Yj) = —G/\ri — rj\ is the Newtonian binary 
potential in c? = 3 dimensions. Starting from the A^-body Fokker-Planck equation and using a 
mean- field approximation [2011^; we can derive the nonlocal Smoluchowski equation 

i(TVp + pV$) 



where the potential $(r, t) is produced by the density p(r, t) according to the Poisson equation 
(3) A<l> = AnGp. 

The self-gravitating Brownian gas has a rigorous canonical thermodynamic structure |,7j. It 
can be seen therefore as the canonical counterpart of the usual A^-stars problem governed by 
Newton's equations and possessing a microcanonical structure. 
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On the other hand, in a biological context, a model of chemotactic aggregation exhibiting 
blow-up phenomena is provided by the system of equations 

(4) ^ = DAp-xV{pVc), 

(5) Ac = -\p, 

where p is the concentration of the biological population (amoebae), c is the concentration of the 
substance secreted (acrasin) and x measures the strength of the chemotactic drift. Equations 
dU and (0) can be viewed as a simplification of the Keller-Segel model Clearly, the two 
models are isomorphic provided that we make the identification $ ^ — ^^c, T ^^^^ and 
^ ^ ^x^' ■'■'^^^o'iucing the mass M = J pd^r of the system and the radius R of the domain, 
we can show that the problem depends on the single dimensionless parameter r] = (3GM/R [7:. 
Therefore, a large value of r] corresponds to a small temperature T or a large mass M. 

2.2 The Smoluchowski-Poisson system 

In the following, we shall use the gravitational terminology but we stress that our results are 
equally valid for the biological problem due to the above analogy. By rescaling the physical 
parameters adequately, it is possible to take M = R = G = $, = 1 without restriction. Then, 
the equations of the problem are 

(6) ^ = V(TVp + pV<l>), 

(7) A$ = 47rp, 

with proper boundary conditions in order to impose a vanishing particle flux on the surface of 
the confining sphere. If we restrict ourselves to spherically symmetric solutions, the boundary 
conditions are 

9$ do 

(8) _(0,t) = 0, $(l,t) = -l, T^(i,t)+p(i,t) =0. 

With this rescaling, the equations only depend on the temperature T, or equivalently on the 
parameter rj = 1/T. As mentioned previously, a small temperature is equivalent (through a 
rescaling) to a large mass. 

Integrating Eq. ((7j) once, we can rewrite the Smoluchowski-Poisson system in the form of a 
single integrodifferential equation 



dt r"^ dr \ \ dr 

The Smoluchowski-Poisson system is also equivalent to a single differential equation 
(10) ^ . T f ^ - 2^ V \m^^' 



dt \ dr"^ r dr J dr 

for the quantity 

(11) M(r,t)= / p{r',t)A7fr'^ dr', 



which represents the mass contained within the sphere of radius r. The appropriate boundary 
conditions are 

(12) M(0,t) = 0, M(l,t) = l. 
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3 Asymptotic estimate of the collapse time 
3.1 Systematic expansion 

We know that the Smoluchowski-Poisson system blows up for T < Tc, where Tc is a critical 
temperature below which there is no equilibrium state [7j. We shall place ourselves close to the 
critical point and expand the mass profile as 

(13) M{r,t) = M,{r) + eMi{r,t) + e^M2{r,t) + 

where Mc(r) is the equilibrium profile at Tc and e is a small parameter defined as 

rj-i rj-i ^ 1/2 



(14) .^^^J «1. 

We also rescale the time according to 

(15) t = r/e. 

If we substitute the expansion Eq. (fT^ in Eq. (fTUI) and equal terms of same order, we get to 
order 0: 



(16) TJ ——--—-] +——- = 0, 



to order 1: 



" \ dr'^ r dr J r"^ dr ^ ^ ' 

and to order 2: 

dMi ^fd'^M2 2dM2 d^M^ 2dMA l/^^^Ma ,,dMi ,,dM^ 



The boundary conditions are 

(19) M„>o(0,r) = 0, M;>o(0,r)=0, 

(20) Me(l,r) = l, M„>i(l,r) = 0. 
3.2 Order 0: the equilibrium state Mc{^) 

At equilibrium, the condition of hydrostatic balance TVp + pV$ = combined with the Gauss 
theorem d^/dr = M[r)/r^ leads to the relation 



(21) M(r) = -Tr 

dr 

Since M'(r) = Anpr'^, we obtain the differential equation 



T =0.39721137... 




In (a) 



Figure 1: We plot r](a) = T^^(a) where a = {Anpo/TY^'^ parameterizes the series of equihbria. 
For a > ac, unstable modes with A > arise, which are solution of Eq. ()33|). The minimum 
temperature Tc corresponds precisely at the point of marginal stability A = 0. 



^o0.5 



T =0.39721137. 

c 



0.2 0.4 0.6 0.8 1 

Figure 2: The integrated mass density is plotted at the critical point. 
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Introducing the function ^p through the relation 

(23) p = poe-^, 

where po is the central density, and using the normalized distance ^ = ar where a = {AupQ/Ty^'^, 
we can rewrite the relation Eq. (PT|) in the form 

(24) M(0=T^^'(0. 

a 



Furthermore, according to Eq. (j22|). the function ip is solution of the Emden equation 



(26) V(0) = ^'{0) = 0. 

The Taylor expansion of the Emden function near the origin is 

(27) ^l: = ^^^ + ... 
^ ^ ^6^ 120^ 

Using these results, we can check that the Emden equation (j25p is equivalent to Eq. as it 

should. For a given temperature T > T^, the Emden equation has to be solved from ,^ = to 
^ = a such that 

(28) ai;'{a) = 1. 

For the profile we are considering, at the critical temperature Tc, we have to stop the integration 
at Oc such that T{ac) is minimum. The condition [1/T]'(ac) = is equivalent to 

It is found numerically that ac = 8.9931149... and = 0.39721137.... 

3.3 Order 1: the function F{^) 

Since Eq. (|17|) is linear, we look for solutions of the form 

(30) M^=A{T)F{0■ 
The function F{^) satisfies the differential equation 

d^F 2dF 1 d . 

which is equivalent to 

, , d fe^dF\ F 
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where we have used the Emden equation Eq. (j^ . Equation is precisely the equation 
found at the critical point Tc by analyzing the linear stability of isothermal spheres 7]. Indeed, 
considering a perturbation SM ~ e'^* around a stationary solution M(r) of Eq. ()10|) . we get 

(33) t(^^-»+1A,m«)^A.M. 

If A < 0, the stationary solution is stable and if A > 0, the stationary solution is unstable. 
At the critical point T^. where A = (marginal stability), we recover Eq. (jl7|) with Mi = 6M. 
Equation (jH^ has to be solved with the boundary conditions 



(34) F(0) = F'(0) = 0. 

We find that -F(0 can be simply expressed in terms of ipiO as P^ : 

(35) F(0 = e'e-^ - eV- 

The condition F{ac) = determines the critical temperature Tc. Using Eq. ()35|). this condition 
is equivalent to 

Comparing this result with Eq. ()29|1 . we find that the criterion of marginal stability (A = 0) 
corresponds to the point ac in the series of equilibria T = T{a) where the temperature is 
minimum ([1/T]'(aj = 0) [13 HI. 

3.4 Order 2: the function A{t) 

Using the previous results, Eq. (fTH|) can be rewritten 



(37) F(Oi(r) = Teal 



d'M2 2dM. 1 d , 



We now multiply Eq. by a function x(^) and integrate between and ac- We choose the 
function x (see below) so as to eliminate the terms where M2 appears. We are thus left with 

(38) A = KA^ + B, 

where 

For future convenience, we note that 

(40) F'{o = ee-H2-m- 

By construction, the function x satisfies the integral relation 



(41) 



d'M2 2dMo 1 5 



X(0 di = 0. 
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Integrating by parts and using M2(0) = M2(0) = M2{a) = 0, we obtain 

(42) x(«c)M^(«„ r) + ^"10 + Q - ^'^ ^ + ^(e^' - l)x}M2(e, r) rf^ = 0. 
Thus, we impose that x is a solution of the differential equation 

(43) ll+f|-#'')$ + l(«''-l)x- = 0. 



dp Vf /if C- 
with the boundary condition 

(44) X(0) = 0. 

It turns out that the solution of this equation can be expressed in terms of the Emden function 
as 

(45) ;^(e) = lF(Oe^ = e - ^'e^ 

with x'(0) = 2/3. The function xiO is represented in Fig. 3 along with the function F{C,). 
Using Eq. (jHUj) . we readily check that x(q^c) = 0. Therefore, the condition Eq. PT|) is satisfied. 
Accordingly, A{t) is determined by Eq. (j38|) . where K and B have to be evaluated using 

Eq. (jnni)- 




8 



3.5 The blow up time 

The solution of Eq. (|HHj) with A{t) — > — oo as r ^ is 
Returning to the original time variable, we get 

t{i-T/T,y/\BKy/^ . 

Then, we find that the blow up occurs for 

(48) t,oii = Uiv - ^c)-^^^ 
with 

(49) '--(w)"- 

We find numerically that K = 62.56038... and B = 0.4716274.... We thus get 

(50) t^u =t,{r]- 7],)-^l\ U = 0.91767702... 

This asymptotic result is compared in Fig. 4 with the exact law tcou{T) obtained by solving 
numerically the Smoluchowski-Poisson system. 



(46) 





1/Z 


1 


tan 










(T-T) 



Figure 4: We plot tcou x (T^ — TY^'^/Tc computed numerically as a function of (T^ — T)^/^, 
which should converge to = 0.91797702.... A quadratic fit in the region (7;-T)i/2 < 0.05 
retrieves the first four digits of t^,. 
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4 Relaxation time estimate for T > Tc 



We now assume T > so that an equilibrium state exists. We look for a time dependent 
solution for M{r,t) which converges exponentially to the equilibrium profile 



(51) M(r, t) = Mrir) + FT{r)e 



't/r 



The purpose of this section is to compute the leading asymptotic form of the relaxation time 
r(T) +00 as T goes to T^. After inserting this ansatz into the equation of motion for M(r, t), 
we obtain 

(52) + (^0' - ^) + «'e-<^FT = -{tT)-^Ft, 
where 

(53) 0(r) = ip{ar). 

We have used r as a coordinate (instead of ^) so as to remain within a fixed interval of space 
< r < 1 as we make the perturbative expansion in a (see below). When T — T^, we have 
a ^ ac, T ^ +00, and -FT(r) — > F(Q!cr), where F{^) is given by Eq. (jSH). We now expand 
Eq. in power of e = by introducing /(^) and /i such that 

(54) FT{r) = F{acr) + ef\a,r) + 0(5^). 
and 

(55) {Tal%)-^ = ^ie + 0{e^). 

The relation between T — Tc and e will be given later. Collecting terms proportional to e, and 
setting ^ = acX, we obtain 

(56) /" + (^^' - ^) /' + /e"'^ - ^FF' = -/iF, 

subject to the boundary condition /(O) = /(etc) = 0. This eigenvalue problem selects a unique 
value for /j,, hence for r. We were not able to solve this problem analytically. However, we 
remark that this equation is similar to Eq. (jHTj) . without the term leading to B. We can 
therefore use the same trick as in Sec. 13.41 After multiplying Eq. ()56|) by the function x defined 
in the previous section, we end up with 

(57) /i = ,f , : \ , ^ = ^ = 0.17202792... 

The direct numerical solution of Eq. leads of course to the same value. The function f{^) 
is plotted in Fig. 3. 

We now find a simple relation between T — and e close to ^31- Using the condition 

aip'^a) = 1/T and the fact that (^^')|^ = O5 "^^ obtain 

(58) ^ - ^ = ^^^^^mi + 0((«c - af), 
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which up to leading order in e can be written in the form 
(59) e ~ 



2(T - Te) 



Te(l - 2T,) 

Finally, we obtain the following diverging behavior for r: 



1 / 2Tc(T — Tc) 2 



(60) 

Numerically, we get 

(61) r ~ ct{T - T,r'/' ~ c,(r/, - r^y'/' 
with 

(62) ct = 0.036563056..., = 0.09204937... 

Note that a heuristic argument leading to r ~ (T — T^)^^/^ was given in (7], although the 
constant ct could not be calculated. 

5 The decay of the density profile at T = 

Strictly at T = Tc, r = +oo, and one expects a slow convergence to the equilibrium profile 
MJr). Let us write 



(63) M(r,t) = M^{r) - A{t)Mi{r) - B{t)M2{r) + ... 

We implicitly assume that B(t) ^ A{t) as t ^ +oo. Moreover, the presence of the minus signs 
in Eq. (j63j) anticipates the fact that starting from a fiat density the central density increases. 
In |7], it was argued based on a heuristic argument that A{t) ~ t~^. It is the purpose of this 
section to justify this statement as well as to give more quantitative estimates. 

Let us postulate for the moment the natural choice B{t) ~ A(t)^ which will be justified 
hereafter. Inserting again this ansatz in the dynamical equation for M{r,t), and collecting 
terms of order A{t) and A{ty, we immediately find that the first equation leads to 



/2T 

(64) Mi(r) = J^F{a,r), 

V etc 

where the constant has been set for later convenience, as multiplying Mi and dividing A by the 
same constant leaves the expansion for M unchanged. Setting M2(r) = g{acr), we find that g 
satisfies 

(65) g" + (^V' -f)9' + 9^-^ - ^^FF' = '^F, 
with 
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which must be a constant independent of time. Thanks to the proper choice of the constant 
in Eq. we find that Eq. ()65p and Eq. (jSUj) are identical. The only physical solution with 
g{0) = g{ac) = corresponds to the eigenvalue A = Finally, we find 

(67) A(t)~t-V-^y^a-l 

Once the signs have been fixed in Eq. (jUHj) . we indeed find that the small time correction to 
the central density is necessarily negative, as Mi, M2 (or g) and A are found to be positive. 
Therefore, the equilibrium profile Mc(r) is approached from below, which is natural since an 
excess of mass would provoke gravitational collapse. 

Let us illustrate this result by calculating the correction to the central density. Using 
Eq. (IHH) and Eq. (jHIj), and the fact that Mc(r) ~ ^por^ and F(^) ~ |,^^, we find that for large 
time the central density converges to po from below in a universal manner 

(68) po-p(r = 0,t) = ^ + 0(r2), 
with 

(69) Cp = (nfi)-^ = 1.8503385... 
This result is illustrated quantitatively in Fig. 5. 




t 

Figure 5: At T = T^, and starting from a uniform mass density, we have computed numerically 
Po — pij" = 0,)!:), where po is the equihbrium central density. This is compared with the exact 
and universal large time estimate Cp/t, with Cp = (vrp)"^ = 1.8503385... (dashed line). 

Finally, we discuss the a priori unjustified choice B(t) ~ A'^(t). The case B(t) <^ A'^it) 
is clearly impossible as it amounts to take g = in Eq. (j65|) . leading to an equation which is 
never satisfied by F, whatever the choice for A{t). The opposite case B{t) ^ A'^{t) leads to an 
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equation similar to Eq. (j65|) but where the term proportional to FF' is absent. For this equation, 
it is clear that if a solution exists for a specific Aq, then X/XqQq is also solution associated 
to A. This implies that the only solution satisfying the boundary condition g{^) = g{ac) = 
is actually the null function associated to A = 0, which is not permitted for T > T^. Hence 
B{t) ~ A'^it) is the only possibility leading to a physically sensible solution. 

6 The pulsation period of isothermal spheres 

In Sec. EJ we have determined the relaxation time of self-gravitating Brownian particles de- 
scribed by the Smoluchowski-Poisson system close to the critical point Tc. Writing the pertur- 
bation as 6p ~ e'^* where A = — 1/r, the eigenvalue equation for A can be written in dimensional 
form as 

^ f ^ d<l\ Gq A^ 
(70) IT + ^ = .^.rr^o Q, 



dr dr J Tr"^ AirpTr'^ 

where g(r) = 6M{r). This equation is the same as Eq. (j52p . Using the result of Sec. 01 and 
returning to dimensional parameters, the inverse relaxation time is given by 

(71) X = T-iv.-vY^'^^ 
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The sign — corresponds to a < ac and leads to exponentially damped pertubations (stable). 
The sign + corresponds to a > ac and leads to exponentially growing pertubations (unstable). 

We now consider the case of gaseous self-gravitating systems (stars) described by the Euler- 
Jeans equations: 

(72) 1^ + V(pv) = 0, 



(9v 1 
(73) _ + (vV)v = --Vp- V<l>, 



(74) A$ = 4nGp. 

We use an isothermal equation of state p = pT and we assume that the system is confined 
within a region of size R. These equations provide the usual starting point for the analysis of 
the gravitational stability of gaseous systems. The case of an infinite homogeneous medium 
was first investigated by Jeans in his classical study. The case of a finite inhomogeneous 
medium was considered by Chavanis In that case, the equation of pulsation can be written 
in the form 

d f 1 dq\ Gq X' 



dr \47rpr^ dr J Tr^ AnpTr'^ 

where, as before, 6p ~ e'^*. Comparing with Eq. (fTOj) . we see that the results of Sec. HI remain 
valid provided that A^ is replaced by A^. Therefore, Eq. ()71|) translates into 

(76) ^' = ^^(^^-^)'''^- 
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For a > Oc, the eigenvalue A = ±vA^ can be positive implying instability. For a < Oc, the 

eigenvalue A = ±iv— A^ is purely imaginary so that the time evolution of the perturbation has 
an oscillatory nature. Close to the critical point, the pulsation period is given by 



In this paper, we have obtained the asymptotic expressions of the blow-up time and relax- 
ation time of self-gravitating Brownian particles and biological populations close to the critical 
point of collapse in d = 3. An excellent agreement is obtained by solving numerically the 
Smoluchowski-Poisson system. This study confirms and improves the results obtained in [7] 
on the basis of heuristic arguments. A possible extension of this work would be to consider 
the case of polytropic distributions arising in case of anomalous diffusion PU] • More generally, 
we could consider the generalized Smoluchowski equation (81) proposed in [T7] which is valid 
for an arbitrary equation of state. This equation can provide (among other applications) a 
generalized chemotactic model valid when the diffusion coefficient or the drift term depends on 
the density. These extensions will be considered in future works. 
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